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We calculate the static polarizability of multilayer graphene and study the effect of stacking 
arrangement, carrier density, and onsite energy difference on graphene screening properties. At low 
densities, the energy spectrum of multilayer graphene is described by a set of chiral two-dimensional 
electron systems and the associated chiral nature determines the screening properties of multilayer 
graphene showing very different behavior depending on whether the chirality index is even or odd. 
As density increases, the energy spectrum follows that of the monolayer graphene and thus the 
polarizability approaches that of monolayer graphene. The qualitative dependence of graphene 
polarizability on chirality and layering indicates the possibility of distinct graphene quantum phases 
as a function of the chirality index. 



Since the discovery of graphene in 2004, experimen- 
tal fabrication of gated graphene, it has been a sub- 
ject of great interest [l|. Recently, multilayer graphene 
has attracted considerable attention for its fundamen- 
tal properties and chiral electronic structure, and for fu- 
ture applications due to the presence of additional layer 
degree of freedom. Interestingly, the electronic struc- 
ture of multilayer graphene is sensitive to its stacking 
sequences [2I, [sj, indicating the possibility of engineering 
electronic properties by tuning the stacking arrangement. 
Optical properties of multilayer graphene have been stud- 
ied experimentally [4] and theoretically [5], showing that 
characteristic peak positions in the optical conductiv- 
ity depend on stacking sequences. Transport properties 
of multilayer graphene depend sensitively on the num- 
ber of layers and stacking sequences in the system [6]. 
Quantum Hall effects in trilayer graphene show different 
quantization rules compared with monolayer and bilayer 
graphene [7] . Electron-electron interactions also play an 
important role in multilayer graphene, especially in pe- 
riodic ABC stacking compared with other stacking se- 
quences due to the appearance of relatively flat bands 
near the Fermi energy indicating the possibility of in- 
teraction induced ordered states [8] which could be very 
different from those studied in monolayer and bilayer 
graphene. 

As the number of graphene layers increases, effects of 
screening become more and more important due to the 
enhancement of the electronic density of states (DOS), 
which determines the fundamental properties such as 
transport, optical properties, phonon dispersion, and 
RKKY interaction. The energy band structure of multi- 
layer graphene is very sensitive to its stacking sequence, 
leading to its screening properties depending strongly on 
the stacking arrangements. There has been extensive the- 
oretical activity on Coulomb screening in monolayer [9j 
and bilayer graphene [10], but there is no systematic study 
whatsoever of screening in multilayer graphene. 

In this Letter, we calculate the static polarizability 
of multilayer graphene and systematically study the ef- 
fect of stacking arrangement, carrier density, and onsite 



energy difference on its screening properties. We first 
study the screening properties of a chiral two-dimensional 
electron system (C2DES) with an arbitrary chirality J 
demonstrating that the polarizability and the associ- 
ated screening behavior strongly depend on J due to 
the enhanced (suppressed) backscattering for even (odd) 
J. Next we calculate the screening properties of multi- 
layer graphene, and then compare and analyze the re- 
sults using those of C2DES. At low densities, the energy 
spectrum of multilayer graphene is described by a set 
of C2DESs[3| and screening is determined by its chiral 
nature. As density increases, the energy spectrum even- 
tually follows that of monolayer graphene and thus the 
polarizability approaches that of monolayer graphene in 
the high density limit. We also show that the onsite en- 
ergy difference between the sublattices or layers, which 
could induce a band gap, can enhance or suppress the 
backscattering depending on the chiral nature of multi- 
layer graphene. 

The static polarizability is defined by 
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where g = g^gv is the total degeneracy factor [g^ = g^ = 2 
are spin and valley degeneracy factors, respectively), ex^k 
and fx,k are the eigenenergy and Fermi function for the 
band index A and wavevector /c, respectively, Fx^x'{k, k') 
is the square of the wavefunction overlap between \\^k) 
and I A', k') states, and k' = k + q. 

To understand the screening properties of multilayer 
graphene, we first consider a 2-band pseudospin Hamil- 
tonian which describes 2D chiral quasiparticles[3]. The 
pseudospin Hamiltonian with the chirality index J is of 
the form 
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2 



Fermi velocity and t± is the nearest-neighbor inter- 
layer hopping. The Hamiltonian has an energy spec- 
trum given by Sg^k = st± ^^^yZ^) corresponding 
eigenfunctions are \s,k) = -^(5,6*'^'^^), where (j)^ = 

tdin~^ (ky/kx) and 5 = ±1 for positive (negative) energy 
states, respectively. Note that for C2DESs, Fg^s'ik, k') = 
^ [1 + ss^ cos J(0fc — ^fc/)], thus the backscattering with 
0fc — = ±7r between the same bands s = s' is en- 
hanced (suppressed) for even (odd) J due to the chiral 
nature of the electronic structure. This has important 
consequences as scattering properties of odd/even J turn 
out to be qualitatively different. 

To investigate the polarizability we rewrite Eq. ^ as 
n(g) = nintra(^) + ninter(^), whcrc the iutrabaud polariz- 
ability Ilintra IS defined by the intraband transition terms 
A = AMn Eq. ([1]), and the interband polarizability Ilinter 
with A ^ A' in Eq. ^ which is induced by the virtual 
interband transitions. 

For the undoped intrinsic case (with the carrier den- 
sity n being zero), we find Ilintra (<?) = and ninter(<?) = 
N{q)I{J), where N{q) = ^^j^lfj^n^y and /(J) is a 
constant depending only on the chirality index J, given 

by 
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where x' = y 1 + 2xcos0 + x^, and cos^ = (x + 
cos(j))/x'. For example, we have /(I) = 7r/8 and /(2) = 
log 4. Note that the screened Coulomb potential is given 
by V'sc(^) = 27re^/ (g + gs(g'))], where is the back- 
ground dielectric constant and qs{q) = {27re'^ / hz)Il{q) is 
the screening wavevector. For the intrinsic case, ^s(^) ^ 
g'^"^, and thus in the long wavelength limit, the screened 
Coulomb potential behaves as V^c(^) oc g^"^ and for 
J > 3 it vanishes as g 0, whereas Vsc{q) oc 1/q as 
g ^ 00 for all J. Thus, the long wavelength Coulomb 
interaction is completely screened for J > 3. 

Next we consider the polarizability for the doped ex- 
trinsic case (i.e. n 7^ 0). At = 0, the polarizabil- 
ity becomes the DOS at Fermi energy sp, i-e. 11(0) = 

7V(£f) = No where Nie) = ^^^^ (klA^)*"' is the 
DOS at energy e for C2DES. Note that in the intrin- 
sic limit (n 0, or equivalently ep 0), Nq diverges 

as ^ or n^~i for J > 3, which indicates enhanced 
screening at low densities. In Fig. [TJa) the Thomas- 
Fermi (TF) screening wavevector, q^F = Qs{q = ^) = 
{27Te'^ / hi)Il{q = 0), is shown as a function of density, in 
which qTF ^ oc 2" diverges as n ^ for J > 3. For 
comparison, we also show the TF screening wavevector 
for multilayer graphene in Fig.fljb). At low densities, the 
energy spectrum of multilayer graphene is described by 
a set of C2DESs[3], and the C2DES with the largest chi- 
rality dominates the DOS and the screening wavevector. 
Note that the periodic ABC-stacked J-layer graphene, 
whose energy spectrum is described by J C2DES at low 
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FIG. 1. (a) The Thomas- Fermi wavevector, qtf oc A^o = 
Il{q = 0), as a function of density n with k = 2.5, where 
a — 0.246 nm is the lattice constant of graphene. (b) 
Same as (a) for multilayer graphene with several stacking se- 
quences, where MLG and BLG represent monolayer and bi- 
layer graphene, respectively, (c) and (d) show the total static 
polarizability n(^) and its first derivative dll{q)/dq^ respec- 
tively, as a function of wavevector for C2DES with chiralities 
J — 1,2, 3, 4, 5. Note that the first derivative at 2/cf is con- 
tinuous (discontinuous) for odd (even) J. (e) and (f) show 
the intraband and interband polarizabilities, respectively. 



densities, shows divergent screening behavior for J > 3. 

Figure [T]^c) shows the total static polarizability as a 
function of wavevector for C2DES. We note that the nor- 
malized polarizabilities of Fig.[T]^c) vary in orders of mag- 
nitude depending on J if they are expressed in absolute 
units, as indicated in Fig. [IJa). In Fig. [Hd) the first 
derivative of the total polarizability with respect to the 
wavevector, dll{q)/dq, is shown for each J. In Figs. [TJe) 
and (f) we show the intraband and interband polarizabil- 
ities for several chirality indices J = 1,2,3,4,5, respec- 
tively. 

For small the polarizability becomes 
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As q increases, Ilintra (Hinter) decreases (increases) 
quadratically in the leading order and these quadratic 
leading order terms exactly cancel out in the total 
polarizability. Thus, the total polarizability increases 
slowly as except for J = 1, where two terms ex- 
actly cancel out up to = 2/cf giving a constant 
2D polarizability for q < 2/cf. For large we have 

Hintra = iVo [{J /2){q/ k^)'^ ^ 0{q-^-^)] and U.^ter = 

No [l{J){q/kF)'^-^ + 0((7"^)]. Thus for large wavevec- 
tors q ^ 2/cf, the polarizability from interband transition 
dominates, and n(g) ^ ninter(^). 

Since many significant physical properties (e.g. trans- 
port, Kohn anomaly, and quantum criticality) are in- 
duced by the behavior of the polarizability dit q = 2/cf, 
we look further into Il{q) near q = 2A:f. At q = 2/cf, the 
calculated polarizability with even J shows a cusp while 
that with odd J varies continuously with q appearing as 
an infiection point in Il{q). As q approaches 2A:f from 
above we find 11 (g) oc (g^ — 4A:p)^/^ for odd number of 
J and Il{q) cx (g^ — for even number of J. As 

a consequence, the static polarizability has a discontin- 
uous (continuous) first derivative for even (odd) J [see 
Fig. [11(d)]. For even Js the first derivative is discon- 
tinuous at g = 2A:f showing dll{q)/dq cx \l ^Jq^ — 4/c|, 
while for odd Js, the discontinuity appears in the second 
derivative as (Pll{q)/dq'^ oc — 4/c|. The origin of 

these singular features is closely related to 2/cf backscat- 
tering in the system. When the contribution of the 2A:f 
scattering to the polarizability is enhanced due to the chi- 
ral property of even J (or the overlap factor in Eq. ([T])), 
the polarizability has a singular feature, whereas when 
the backscattering is suppressed for odd J, no singular 
features occur. Note that all these singular properties 
arise only from the intraband polarizability. 

So far we have studied the polarizability for a simple 
two band model of 2D chiral quasiparticles. Now we con- 
sider the polarizability for a full band continuum model 
of multilayer graphene near the hexagonal corners of the 
Brillouin zone, called the K and K' points, taking into 
account only nearest-neighbor intralayer hopping t and 
the nearest-neighbor interlayer hopping t_\_. In this pa- 
per t = 3 eV and t± = 0.3 eV will be used in numerical 
calculations. 

Figure [2fa) shows the polarizability of zero-gap bilayer 
graphene as a function of wavevector for several carrier 
densities with the energy band structure as an inset. For 
small densities, the polarizability of the full four band 
model resembles that of the J = 2 C2DES (top black 
line). As the density increases, the polarizability eventu- 
ally approaches that of monolayer graphene because at 
high densities the interlayer coupling becomes negligible 
and the energy band structure behaves like a collection 
of monolayer graphene sheets. The main difference be- 
tween the polarizability of monolayer and high density 
bilayer graphene is the existence of a kink at q = 2A:f, 
which arises from the enhanced backscattering and does 
not vanish completely in bilayer graphene even at very 
high densities, i.e. the system maintains its memory even 




FIG. 2. (a) Static polarizability of bilayer graphene with the 
zero energy gap as a function of wavevector for several carrier 
densities in units of 10^^ cm~^. The inset of (a) shows the 
energy band structure of bilayer graphene. In (b) we include 
the onsite energy difference U between the layers to open 
the energy gap. The polarizability is calculated for a carrier 
density n = 10^^ cm~^ and for different values of U. In both 
figures the black solid line corresponds to the result of C2DES 
with J — 2. The inset of (b) shows the static polarizability of 
monolayer graphene with energy differences between the two 
sublattices of [/ = 0, 0.1, 0.2, 0.3 eV. 



at very high carrier density where the energy dispersion 
has become the same as that in monolayer graphene. 

Figure Efb) shows the static polarizability of bilayer 
graphene with an onsite energy difference between the 
layers for a fixed density n = 10^^ cm~^. As the on- 
site energy difference increases, the energy band gap also 
increases, reducing the interband contribution and the 
overall magnitude of the polarizability. Note that in the 
presence of the onsite energy difference, the backscat- 
tering is suppressed and thus the cusp structure at 2/cf 
becomes weakened. In case of monolayer graphene, how- 
ever, the onsite energy difference between the two sub- 
lattices gives rise to the enhanced backscattering in addi- 
tion to the opening of an energy gap. As a consequence, 
the polarizability of monolayer graphene with the onsite 
energy difference develops a cusp structure at 2/cf, as 
shown in the inset of Fig. [2fb). This can be understood 
if we consider the effect of the onsite energy difference 
on the overlap factor in Eq. ([1]), which weakens the en- 
hancement (suppression) of the backscattering for doped 
bilayer (monolayer) graphene. Therefore, depending on 
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FIG. 3. Calculated polarizability as a function of wavevector 
with several densities in units of 10^^ cm~^ for (a) ABC tri- 
layer graphene and (b) ABA trilayer graphene. The black line 
in (a) and (b) indicates the result of C2DES with J = 3 and 
J = 2, respectively. Insets show the energy band structure of 
each stacking. 



the chiral properties of multilayer graphene, the onsite 
energy difference may enhance or weaken the backscat- 
tering and the cusp structure at 2/cf. 

In Fig.[3]the polarizability of trilayer graphene is shown 
with the energy band structure as an inset. Exclud- 
ing consecutive stacking sequences, there are two dis- 
tinct stacking arrangements, ABC and ABA in trilayer 
graphene. Figure [3fa) shows the polarizability of ABC 
trilayer graphene as a function of wavevector. Note that 
for ABC stacking the low-energy spectrum is described 
by J = 3 C2DES. For comparison, we also show the po- 
larizability of J = 3 C2DES (top black line). At low 
densities, the polarizability resembles that of the low- 
energy C2DES with J = 3 for ABC trilayer, while as 
the density increases, it approaches that of monolayer 
graphene. Figure ^h) is the same as Fig. [3l^a) but 
for ABA trilayer graphene. Note that at low densities 
ABA trilayer graphene is described by a direct product of 
J = 1 and J — 2 C2DESs and the Fermi energy crosses 
two energy bands at k^^ and /cf2 corresponding to the 
J = 1 and J = 2 energy bands, respectively. Thus the 
ABA trilayer graphene has two different Fermi wavevec- 
tors which can be calculated from the carrier density 
n = £-kl = ^(/cpi + )• The polarizability of ABA 
trilayer graphene in Fig. [^b) appears to be qualitatively 
different from that of ABC trilayer graphene. One sig- 




FIG. 4. (a) Density oscillations of C2DES with chiralities 
J = 1,2,3,4,5 using RPA for n = 10^^ cm-^ with k = 2.5. 
In (b) Sn(r) is rescaled by where /3 = 2 (/3 = 3) for even 
(odd) J, which clearly shows the decaying pattern at large 
distances with the periodicity of tt/Zcf- 



nificant feature is the existence of a cusp at 2/cf2 instead 
of 2A:f. As the density increases, the cusp becomes weak- 
ened and it shifts because the ratio /cps/^F decreases. 
This peculiar feature can be understood from the DOS of 
ABA trilayer graphene. Since the DOS of J = 2 C2DES 
is larger than that of J = 1 C2DES at small energies (i.e., 
N{e) oc e^~^) many physical properties of ABA stacking 
trilayer graphene are determined by J = 2 C2DES. Thus 
for low densities the polarizability resembles that of the 
J = 2 C2DES and the cusp structure appears at 2/cf2- 

Similar analysis can be applied to arbitrarily stacked 
multilayer graphene and we can easily generalize the re- 
sults of bilayer and trilayer graphene discussed so far. At 
low densities, the energy spectrum of multilayer graphene 
is described by a set of C2DESs and the polarizability 
behaves as that of a C2DES with the largest chirality. 
As the density increases, the polarizability follows that 
of monolayer graphene because the energy spectrum be- 
haves like a collection of monolayer graphene sheets. In 
the presence of the onsite energy difference, the backscat- 
tering may be enhanced or suppressed depending on the 
chiral nature of the system, giving enhanced or sup- 
pressed cusp structure at 2/cf. 

One direct consequence of the polarizability function 
we have calculated is the density fluctuation induced by 
a point charge which is given by 



6n{r) = J ^Mqr) [e 



(5) 



where Jo{x) is the zeroth-order Bessel function and e{q) 
is the static dielectric function, which is given by e{q) = 
1 + ^^|-n(qr) within the random phase approximation 
(RPA). Figure m^a) shows the calculated density oscilla- 
tions of C2DES for n = 10^^ cm-^ with = 2.5. At large 
distance {k-pr ^ 1) we find a density oscillation (Friedel 
oscillation) which has the form sin(2/cFr)/r^ for even J 
while cos(2/cFr)/r^ for odd J, as seen in Fig. SJ^b). This 
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radial decaying pattern arises from the kink structure at 
2/cf of the polarizabihty function. When the discontinu- 
ity appears in the first (second) derivative of the polar- 
izability, (ni{q)/dq {d?Ii{q) / dq^) ^ the radial decay has a 
form of (r~^). The discontinuity of the derivative 
is the direct consequence of the enhanced (suppressed) 
backscattering for even (odd) J. One interesting aspect 



of the polarizabihty of C2DES, which is reflected in these 
Friedel oscillations as well as in the screening properties, 
is that even (odd) J behaves like qualitatively similarly to 
the Lindhard polarization function for 2D (3D) electron 
gas systems. 
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